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A STUDY OF THE NONLINEAR AERODYNAMICS OF BODIES 
IN NONPLANAR MOTION* 

Levis Barry Schiff 
Ames Research Center 

SUMMARY 

Concepts from the theory of functionals are used to develop nonlinear 
formulations of the aerodynamic force and moment systems acting on bodies 
in large-amplitude, arbitrary motions. The analysis, which proceeds for- 
mally once the functional dependence of the aerodynamic reactions upon the 
motion variables is established, ensures the inclusion, within the result- 
ing formulation, of pertinent aerodynamic terms that normally are excluded 
in the classical treatment. Applied to the large-amplitude, slowly varying 
nonplanar motion of a body, the formulation suggests that the aerodynamic 
moment can be compounded of the moments acting on the body in four basic 
motions: steady angle of attack, pitch oscillations, either roll or yaw 

oscillations, and coning motion. Coning, where the nose of the body 
describes a circle around the velocity vector, characterizes the nonplanar 
nature of the general motion. 

With the above motivation, a numerical finite-difference technique is 
developed for computing the inviscid flow field surrounding a body in con- 
ing motion in a supersonic stream. Computations carried out for circular 
cones in coning motion both at low supersonic and hypersonic Mach numbers 
confirm the adequacy of a linear moment formulation at low angles of attack 
At larger angles of attack, however, the forces and moments become non- 
linear functions of the angle of attack. Computational results for the 
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reactions on the circular cone at the higher angles of attack agree -well 
with experimental measurements within the range of variables investigated. 
This indicates that the initial nonlinear behavior of the aerodynamic 
forces and moments is governed primarily by the inviscid flow. 

1. INTRODUCTION 

Linear formulations of the aerodynamic .force and moment systems do not 
properly describe the aerodynamic reactions on flight vehicles in nonplanar 
motions at large resultant angles of attack. As a result, equations of 
vehicle motion incorporating linear aerodynamic formulations have often 
failed to correctly predict the variety of motion such flight vehicles can 
experience. In the past this has been a problem primarily associated with 
the flight of slender bodies of revolution. However, the requirements for 
increased angle-of-attack ranges for proposed high performance aircraft and 
space shuttle vehicles, as well as those envisioned for STOL aircraft, have 
tended to make the deficiencies of the linear formulation a problem of more 
widespread concern in the analysis of vehicle motions. 

In the linear formulation, a reference flight condition is chosen, 
for example, the steady level flight of an airplane, and the deviations of 
the angular orientation and angular velocities of the body are measured 
from the reference state. The aerodynamic reactions are expressed in a 
Taylor series expansion in terms of the deviations and only terms linear 
in the disturbance quantities are retained. The coefficients of the expan- 
sion are called the aerodynamic derivatives and are evaluated at the refer- 
ence state. When applied to a combined pitching and yawing motion, the 
linear formulation can be shown to be equivalent to the vector decomposition 
of the nonplanar motion into two orthogonal planar motions, to the subsequent 
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treatment of each planar motion as if the other we re absent, and to the 
superposition of the results. This approach has had a great deal of suc- 
cess, particularly for the case of airplane motions where the deviations 
from the reference state are small. A rather complete treatment of the 
linear formulation has been presented by Etkin . 1 At large resultant angles 
of attack, however, it is physically clear that the reactions due to 
motion in one plane will be influenced by the presence of the other motion, 
and thus a more precise formulation will be necessary to account for this 
coupling. 

The form that extensions of the linear formulation should take to 
account for the large angular deviations from the reference state has not 
yet been settled. Guided by the fact that the static forces acting on a 
body of revolution lie in the plane of the resultant angle of attack, 
Nicolaides et at . 2 and Ingram , 3 concerned with missile aerodynamics, 
assumed that the form of the nonlinear generalization for a body of revo- 
lution was the same as that of the linear formulation, but that the aero- 
dynamic derivatives were nonlinear functions of the magnitude of the 
resultant angle of attack. If this formulation is applied to the combined 
pitching and yawing motions of a nonspinning axisymmetric body, it predicts 
that the aerodynamic damping in the plane of the resultant angle of attack 
is equal to that acting perpendicularly to the angle-of-attack plane. It 
can be shown by a comparison of the experimental results of Iyengar 1 * with 
those of Schiff and Tobak 5 that this is untrue for such bodies at large 
angles of attack. Murphy 6 proposed an extension of the linear formulation 
which allowed for the possibility of unequal aerodynamic dampings in and 
normal to the angle-of-attack plane. Unlike the previous one, Murphy’s 
formulation is therefore capable of correctly distinguishing between the 
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out-of-plane damping and the out-of-plane classical Magnus forces in the 
case of the nonplanar motion of a spinning body of revolution. 

Another approach has heen developed by Tobak et at. , 7 “ 10 vho used con- 
cepts from nonlinear functional analysis to develop a formulation of the 
aerodynamic force and moment system for an arbitrarily shaped body that 
does not depend on a linearity assumption. This formulation has heen shown 
to be equivalent to that of Murphy for the special case of a body of revo- 
lution , 9 and reduces to the form of the linear formulation for small resul- 
tant angles of attack. The formulation suggests that the aerodynamic reac- 
tions on a body in an arbitrary nonplanar motion can be compounded of the 
forces and moments acting on the body in four characteristic motions, three 
of which are well known. The fourth, coning motion, in which the nose of 
the body describes a circle around the velocity vector, is seen to have 
particular significance since the nonlinear behavior, with increasing angle 
of attack, of the contribution to the total force and moment due to coning 
motion cannot be evaluated from the contributions due to any planar motions . 
Experimental evaluations 5 5 8 of the contribution due to the coning of a body 
of revolution have confirmed this and have shown this contribution to be a 
potential cause of circular limit motions at large resultant angles of 
attack. In addition, it is anticipated that the contribution due to con- 
ing motion will be important in correctly describing the pre- and post- 
stall behavior of aircraft-like bodies at large angles of attack. 

The objectives of the present work are twofold. The first is to 
review and unify the development of the nonlinear formulation proposed by 
Tobak et aX* and to remove from this analysis an unnecessary assumption of 
constant flight speed. The second, and more important, objective is to 
present a numerical method for computing the flow field surrounding a body 
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in coning motion in a supersonic stream. A finite-difference scheme of 
MacCormack , 1 1 developed as a shock-capturing technique for computing com- 
plex, steady, three-dimensional, inviscid flow fields by Kutler and 
Lomax , 12 is extended to the case of coning motion. The capabilities and 
limitations of the method are described. Results of computations for 
slender conical bodies in coning motion at various supersonic Mach numbers 
are presented and are compared with experimental results and with the 
results of other analytical and numerical methods, where applicable. The 
results will be seen to exhibit significant nonlinear behavior with 
increasing resultant angle of attack, and the significance of the non- 
linearities will be discussed. 

The author wishes to acknowledge and thank Professor Samuel McIntosh, 
Jr., and Professor Holt Ashley for their advice and encouragement during 
the course of this work. Grateful acknowledgment is also given to Murray 
Tobak of the NASA Ames Research Center for his valuable advice and helpful 
consultation, and to Susan Schiff for her encouragement and for her help 
in preparing the manuscript. 

Finally, acknowledgment is given to the National Aeronautics and 
Space Administration for support of the research and of the author’s 
graduate study through the Honors Cooperative Program with Stanford 
University. 
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2. NONLINEAR FORMULATION - REVIEW AND EXTENSION 


In a series of papers, 7 * -10 concepts from the theory of functionals 
were used, to develop a nonlinear formulation of the aerodynamic force and 
moment system acting on a body performing motions of interest, the first 
being the planar motion of an airplane at large angle of attack. 7 The 
analysis was extended to the large-amplitude, nonplanar angular motions of 
a body of revolution whose mass center traversed a straight-line path, 8 
and showed that the total moment could be compounded of the contributions 
from four simple motions. Further extensions of the analysis to the free 
flight of a body of revolution 9 and of an arbitrarily shaped body 10 showed 
that even in these more general cases the total moment still could be 
determined from the contributions from the same four simple motions. The 
resulting formulation allows the angular deviations of the body to be 
large, but is valid only for the low angular rates typical of aircraft 
motions. Unfortunately , a uniform notation was not employed throughout 
the series, while the assumptions of the analysis, covered in detail for 
the planar case, were abbreviated in the later works. Additionally, it 
was unnecessarily assumed that the flight speed remained constant over the 
course of the motion considered. In this chapter the development of the 
nonlinear moment system is reviewed for the large-amplitude planar motion 
of an arbitrary body whose flight speed varies, and it is indicated how 
the formulation can be extended to more accurately represent motions of 
higher frequencies. The formulation is then developed for the most general 
case, that of an arbitrarily shaped body in free flight, again removing the 
restriction of constant flight speed. Finally, the resulting formulation 

is specialized to the previously reported cases. 
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2.1 Coordinate Systems 


Three axis systems having a common origin at the body's mass center , 
and a common axis x^ aligned with a longitudinal axis of the body, will 

D 

be used to describe the motion in all cases considered. Some latitude 
exists in the choice of the longitudinal axis, and this freedom can be 
used to simplify the description of a particular motion. For example, 
when describing the motions of an airplane-like body, the x^ axis is 
usually chosen to be initially coincident with the direction of steady 
flight. Alternatively, when describing the flight of a body of revolution, 
the x.g axis is most often, but not necessarily, chosen as the axis of 
axial symmetry. 

Axes x_, y_, z^ are fixed to the moving body (Fig. l). The compo- 

.D IJ .D 

nents of the flight velocity vector of magnitude V, resolved along the 
x c , y B , z B hody-fixed axes, are Ug, v^, respectively. Thus 


V = 



+ v. 


B 


( 2 . 1 ) 


The resultant angle of attack a is defined by the flight velocity vector 
and the axis. The plane formed by the y^, axes is the crossflow 

plane, illustrated in Fig. l(a). The projection of a unit vector in the 
flight velocity direction onto the crossflow plane is a vector with magni- 
tude <5, and will be called the dimensionless crossflow velocity vector. 
Reference to Fig. l(b) gives 



(2.2a) 
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Y =s — =5 cos a (2.2b) 

e = — = tan a (2.2c) 

Y 

The components of the angular velocity vector of the hody relative to 
inertial space, resolved along the x B , y-g, z B "body axes, are p B> q B , 
r , respectively. 

-LJ 

A second axis system, x , y, z, is chosen to he nonrolling with 
respect to inertial space. Specifically, the component of the angular 
velocity vector of the Xg, y, z axes measured with respect to inertial 
space, resolved along the x^ axis, is zero, while the components 
resolved along the y, z axes are q, ?, respectively. The nonrolling 
axis system has heen used extensively in the field of missile aerodynamics, 
since its use, together with the assumption of small angular deviations, 
leads to closed-form solutions of the equation of vehicle motion. The 
angle $ through which the body axes have rolled at any time t can be 
defined relative to the nonrolling axis system as 



(2.3) 


The angular inclination X of the crossflow velocity vector 6 is mea- 
sured relative to the nonrolling axis system, while ^ is the angular 
inclination of the body axes from the crossflow velocity vector. With 
the aid of Fig. 1(a)* the body roll rate is seen to be the sum 
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P B = X + ^ 


(2.4) 


The components of the angular velocity of the "body relative to inertial 
space resolved in the nonrolling axis system, pg, q, r, are related to 
those resolved in the hody axis system, p B , q^, r^, through 

q. + ir = e " L ^( < l B + (2.5) 

Finally > axes x^ 5 y, z will lie called the aerodynamic axes. Axis 
z lies in the crossflow plane and is coincident with the crossflow veloc- 
ity vector; axis y lies in the crossflow plane normal to the direction 
of 6. The components of the angular velocity of the lody resolved in 
the aerodynamic axis system, q, r, are related to those resolved in 

the lody axes through 

q + ir = e 1 '^(q B + ir B ) (2.6) 

In accordance vith Ref. 10, w_/V will he called the angle-of-attack 

-D 

parameter a, and v B /V will be called the angle-of-sideslip parameter 
g. They are related to the standard NASA definitions of angle of attack 
a and angle of sideslip g through 


tan a = — = — (2.7a) 

U B Y 

V B * , 

sin 0 = “ = g (2.7b) 
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and to 6 and ip through 


a + i& = 6e 1I,J (2.8) 

The components of the aerodynamic force coefficient resolved along 
the y^, 2-n “body axes are CL., C__ s CL, respectively * vhile the cor- 

responding components of the aerodynamic moment coefficient (acting about 
the mass center) are , C^, C^, respectively. Analogously, the compo- 
nents of the force and moment coefficients resolved along the aerodynamic 

axes X-q, y, z are CL, C v , CL and C , C , C , respectively. The compo- 
d A i Zi £ in n 

nents of the aerodynamic moment coefficient resolved in the aerodynamic 
axis system are related to those resolved in the body-fixed axes through 

C 4 = C £ (2.9a) 

C + iC = e U (C rr + iC ) (2.9b) 

m n m n 

The corresponding relations between the components of the aerodynamic 
force coefficient are obtained by replacing £, m, n by X, Y, Z, 
respectively, in Eq. (2.9). 

To completely describe the state of a six-degree-of- freedom free- 
flight motion, it is necessary to specify the velocity and angular velocity 
vectors of the body. These may be expressed in terms of their scalar com- 
ponents resolved in the body- fixed axes u_, v^, p_, q_, , r D , or 

B B ]3 B B n 

equivalently by a, 0, V, p^, q^, r^. Analogously, in the aerodynamic 

axis system, the motion is specified by the scalar variables d , t|>, V, p^, 

• • » 

q, r, or by d, iJj, V, X, q, r since p^ is related to X and ip through 
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Eg* (2.U)* In those motions -where the mass center traverses a straight- 
line path (i.e*, with no lateral plunging), it is easy to show that, in 
the aerodynamic axes , 

<1 = c a - S (2.10a) 

Y 

r = 6$ ( 2 . 1015 ) 

X = yip (2.10c) 

♦ 

where <j> is the coning rate of the resultant angle-of-attack plane around 
the velocity vector* When the plunging of the mass center is eliminated, 
two of the motion variables can be expressed in terms of the regaining 
four* In the aerodynamic axes* a nonplunging motion can thus be described 

by the variables 6, V, (j> since knowledge of 5 and \p implies knowl- 

• • 

edge of 5 and The relations between valuables in the body- fixed axes 
for the case of zero plunging, corresponding to Eg* (2*10), are 


<1 B = "7 + (2.11a) 
r B = " y ^ “ “V (2.11b) 
P B = y 4> + ^ (2.11c) 


Thus it can be shown that, in the body- fixed axes, a nonplunging motion 
can be described by the variables a, g, V, 4>. 

When, in addition to eliminating the plunging of the mass center, the 
angle of attack o, the bank angle the coning rate, and the flight 
speed are all held fixed, the nose of the vehicle describes a circle 
around the velocity vector. This motion (a = const, i|j = const, V = const, 
4> = const) will be called steady coning motion. In the case of coning 
motion, one sees, with the aid of Eq. (2.1l), that 
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P B = Y* 

q B = 6* 
r B = oil 


(2.12a) 

(2.12b) 

(2.12c) 


2.2 Development for Planar Motion 
To illustrate the ideas behind the development of a nonlinear aero- 
dynamic force and moment formulation * we consider, for simplicity, the 
large-amplitude planar oscillations of an aircraft as shown in Fig. 2. 
Assume that prior to time zero the aircraft is in steady level flight. 

At time zero it begins a longitudinal planar motion such that and w^, 
the components of the flight velocity vector resolved along the x^, 2g 
body- fixed axes, respectively, and the angular velocity of the aircraft, 
all vary, while the wings remain level. Thus v^, the component of 
the flight velocity resolved along the y^ axis, and the angular velocity 
components p and r all remain zero throughout the maneuver. The alti- 
tude excursions of the aircraft are assumed to be small enough for the 
atmospheric temperature and density to be considered constant. Further, 
the variation in the total flight speed is assumed to be small enough for 
the effect of Reynolds number variation on the aerodynamic reactions to 
be negligible. Under these conditions the aerodynamic force and moment 
acting on the body at time t after the start of the motion are dependent 
solely on the velocity components u^ and Vg, on the angular velocity 
and on all values taken by these variables over the time interval from 
zero to t. 


2.2.1 Concept of a Functional 

The fact that the aerodynamic reactions on the body at time t are 
dependent not only on the instantaneous values of the motion variables, 
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but also on the past history of the motion, can be expressed mathematically 
by introducing the concept of a functional. 13 Focussing specifically on 
the pitching-moment coefficient (the development for the other force and 
moment components is analogous), one designates the coefficient C m (t) as 
a functional of u^, v^, q^, or alternately as a functional of 5, V, q 
by the use of the square bracket notation introduced by Volterra: 

C m (t) = E[u b (c),w b (0,c 1b ( 5)] = E»[S(0,V(0,q. B (5)] (2,13) 

where 5 is a variable in time running from zero to t . The alternate 

A 

designation is possible since in the case of planar motion, where g, p , 

-D 

r^ are all zero, a and V are related to u^ and w^ through Eqs. (2.l) 
and (2-7) • 

In brief, just as an ordinary function f(x) assigns a number to each 
x for which it is defined, a functional G[y(£)] assigns a number to each 
function y(?) of the set of functions (all of which are defined in some 
interval a < £ Z b) for which the functional is defined. Thus Eq. (2.13) 
may be interpreted as follows: Given any triad of functions a(5), V(0, 

q (O out of the collection of all such triads defined in the interval 

.D 

0 < ^ t, the functional E f assigns a number to 6 (t). 

2.2.2 Nonlinear Indicial Response 

Following Ref. 7, one defines the nonlinear indicial pitching-moment 
response as illustrated in Fig. 3. As shown for the case of a step change 
in V/V^ (where is a constant reference speed), two motions are con- 

sidered. The first begins at time zero, and at time t the motion is 
constrained so that the motion variables a(0* V(c), and q (£) are held 
constant thereafter. The second motion differs only in the step imposed 
at time t. The difference between the pitching-moment coefficients 
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measured at time t for each of the two motions is divided by the magni- 
tude of the step. The limit of the ratio, as the step size approaches 
zero, is defined as the indicial pitching-moment response at the time t 
per unit step in V/V R at time t. The indicial responses to step changes 
in a and q R are defined analogously. In functional notation these are 


AC (t) 

, l im . ~ ( 7 7v ~T = C [S(€),V(g) a q (g);t, T ] <2.lUa) 

A(V/V r H) ^ ' V V * 

AC (t) 

1 i m “fjj — = Cm /v [a(S),V(0,q^(?);t,r] (2.l4b) 

Aa-*-0 a 


lim 

A(q B t/V R hO 


AC (t) 

aTOTO = % ^(U,v(0,<i B U);t,x] 

-B K B 


(2. l4c) 


As indicated by the notation, the indicial response is dependent not only 

on the levels a(x), V(t), and q^i) at vhich the steps occur, but may 

depend on all values taken by them over the time interval zero to t* 

Thus, in general, the indicial responses may themselves be functionals* 

Breaking the time histories of the motion variables a, V, and q_ 

B 

into a series of step changes and summing the incremental responses to 
each of the steps over the interval from zero to t gives an exact inte- 
gral form for C (t): 
m 

£ m (t) = a m (0) + f (S)dx 

rt 

+ J C my [S(c),V(0,q B (c);t,T] (V/V R )dx 

J t 

»V(0»q B (e) ;t ,x] ~ (q B i/V R )dx 


ll+ 


(2.15) 



2.2.3 Approximate Formulation for Slowly Varying Motions 

While exact in principle. Eg. (2,15) cannot generally he evaluated 
since the functionals appearing within the integrals may depend on the 
entire past history of the motion. It is desired to develop an approxi- 
mate, more easily evaluated form of the expression in which the functionals 
are replaced by ordinary functions. To do so, it is first seen from phys- 
ical considerations that the indicial response to a step at time t is 
not equally affected by all past events, but is most strongly influenced 
by events of the most recent past. This is most easily illustrated for 
supersonic flight, where the flow field around the aircraft at the time 
t, and thus the subsequent indicial response, can only be influenced by 
events that have occurred during the time interval i-t a < £ < t. The 
characteristic time t a is approximately the maximum time required for a 
fluid element to travel from the bow shock wave to the tail shock. The 
indicial response to a step at time t is unaffected by events that have 
occurred prior to x-t a , In subsonic flight the vorticity shed by the 
aircraft will influence the local flow field for all later times , but the 
magnitude of the influence decreases as the vorticity is convected far 
downstream. In this case t a is the time required for the shed vorticity 
to have negligible further effect on the flow field. Additional discussion 
concerning the diminishing dependence on the past can be found in Ref. 7. 

The periods of oscillatory aircraft, motions are typically large com- 
pared to the characteristic time t a . For these slowly varying motions, 
the simplest way of accounting for the dependence of the indicial response 
on the past is to assume that the response is dependent only on the levels 
ct(x), V(t), and q^(x) at which the steps occur. As shown in Fig. k for a 
step in a at i, this assumes that the response for motion B in which 
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a, V 9 and q^ are held fixed over the time interval x-t a < £ < r which 
can influence the response is a good approximation of the response to the 
response for the arbitrary motion A. Under this assumption the response 
is no longer a functional, hut is a function of the parameters a(t), 
V(t), and q^fi) and of the variables t and t. In addition, with the 
motion held fixed prior to the step, the response is independent of the 
time at which the step is made and is seen to be a function only of t-x 
rather than of t and t separately , This is expressed in functional 
notation as 


CmJaU)>V(c) s q B U);t,T] ~ C m ^[a(x ) >V(x ) >q. B ( x ) ;t , x] 
a a 

= C m ^(a(x) J V( T ),q B (x);t- T ) ( 2 . 1 6 ) 

a 

Under this assumption, the desired approximate form of Eq. (2.15) is 

- - C ri 

c (t) = C (0) + 1 6 m (ct(x) ,V(t) 9 q.p(x) ;t-x) (a)Ux 

m ml AU ^ n cIt 

J o a 

+ J C rnv (S(x) s V( T ) s CL B (x);t-x) (V/V R )d T 
/•t 

+ J ^m^ (a(x) ,V( T ) ,q B (x) ;t- T ) ^ (q B £/V R )d T 

(2.17) 

When the deficiency functions (which vanish for t-x > t a ) are introduced 
and defined as 


F(a(x ) ,V(x) »0. B (x ) jt-x) = C m ^(a(x) ,V(x),q B (x) ;«>)-C m (a(x) ,V(x) ,q B (x) ;t-x) 

a a 

(2.18a) 

G(S(t) ,V( x) ,q B (x ) ;t-x) = ^^(aCx ) ,V(x ) ,q R (x) ;»)-6 mY (o(x) »V(x) ,q B (x) ;t-x) 

(2.18b) 


H(a(x) ,V(x ) ,q B (x ) ;t-x) 



(2.18c) 


1 6 



where C m ^ (£( t ) ,V(t ) , q^ ( t ) ;«) is the steady-state value of the indicial 
a 

response, Eq. (2,17) can be expressed in the alternate form: 


e .<‘> 


= 6 m («;a(t) ,V(t) ,q B (t) ) - j* F(S(t ) ,V(x ) ><1 b (t ) ;t— t ) ~ (a)dT 

J Q 

-f 

J o 


G(a(T),V(T),q. B (T);t-t) (V/V R )dx 


'I 


H(S(t) ,V(t) ,q B (x) ;t-r) — (q. B Jl/V R )dT 


(2.19) 


since the terms involving the steady-state values form a perfect differen- 
tial. The term ^ m (°S“(t ) ,V(t) ,q. B (t ) ) is the steady pitching-moment coef- 
ficient that would be measured in a steady pullup maneuver in which a, V, 
arid q g are held fixed at S(t), V(t), and q^ft), respectively* Equa- 
tion (2*19) can be further simplified consistent with the assumption of 
slowly varying motions used to approximate the functional indicial 
responses. The functions are expanded in a Taylor series about the point 
a - a(t) s V ~ V(t) 5 q - 0, and only first-order frequency terms are 
retained. The resulting nonlinear formulation is 


C m (t) = C m (“;a(t),V(t),0) + 


q. B (-t ) 


c m (« ; a(t) s v(t),o) 

R q B 


+ C m .(a(t),v(t),0) + hill C m -U(t),V(t), 0 ) 

R “ V R " 


(2.20) 


where 


v r r t 

x = - — \ F(a(t),V(t) 5 0; T )dT 

a J o 

V R * 

=-~l G(a(t),V(t),0;T)dT 

* o 


Pm* “ 
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Each term in Eq* (2, 20) is again identified with a particular motion from 
which it may he evaluated* The first term is the static pitching coef- 
ficient measured with a and V held fixed at their instantaneous values, 
with q^ held fixed at zero* With the reference velocity V chosen to 
he equal to the instantaneous value of the flight velocity V(t), the 
second term is seen to he the rate of change with q^, evaluated at 

D 

q = 0, of the pitching-moment coefficient that would he measured in a 
steady vertical pullup with q^ = const, a and V held fixed. The third 
term is the damping coefficient that would be measured for small oscilla- 
tions in a about the instantaneous value a(t) = const, with V held 
fixed at V(t), q^ = 0, a plunging oscillation normal to the velocity 
direction. The last term is the damping coefficient measured for small 
oscillations in V about the fixed flight speed V ( t ) , with S held 
fixed at a(t), q = 0, i.e. , the damping due to lunging oscillations of 
the mass center along the flight velocity direction. Thus the specifica- 
tion of the moment acting on the aircraft performing an arbitrary planar 
motion requires knowledge of the moment acting on the body in four 
characteristic motions. 

The number of characteristic motions required can be reduced in 
those cases where the lateral plunging of the mass center is small. To 
demonstrate this we rewrite Eq. (2*20) as 


6 m (t) = e m (“^(t),v(t),o)+ (“‘;S(t),v(t),o)+ Y e m .{S(t) f v(t),o)j 




*B 


VftU - 1 ( Yl-ofa )“S(t ) ) i 

+ - 1 C m .(S(t),V{t),0)+ i G m (»;S(t),V(t),0) 


B 


A B 


( 2 . 21 ) 
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The terms in Eq. (2.21) are identified by comparing them with those 
obtained for the case of nonplunging motion at constant flight speed* 


Under these conditions* V = 0* q_ = (l/y)a, and the last two terms of 
Eq. (2.21) vanish* With the plunging of the mass center eliminated, 
changes in a correspond to angular motion of the body about the fixed 
y_ axis. Thus the term (C m + yC^.) is seen to be the planar pitch 

CL 

damping coefficient that would be measured in a single-degree-of- freedom 
experiment in which the body performs small angular oscillations about a 

mean angle of attack, i.e*, small oscillations in a about a = const, 

* 

with V held fixed and ri = (l/y)a* In a general motion, the contribu- 

Ij 

tions from the last two terms of Eq. (2.2l) are not zero. However, for 
motions in which the plunging is small, and where the flight speed makes 
only small oscillations about a mean speed, the contributions from these 
terms can be justifiably neglected since, in equations of vehicle motion, 
they would appear only as products of (relatively small) damping terms. 

In such cases the total moment acting on the aircraft is due to the con- 
tributions from only two characteristic motions: steady angle of attack 

and damping-in-pitch. Here V(t) need not appear explicitly within the 
notation, it being understood that the characteristic motions will be 
evaluated at a fixed speed equal to the mean value of the flight speed. 


2.2.4 Extensions to Describe More General Motions 


To obtain the approximate integral aerodynamic formulation, Eq. 

(2. IT), from the exact functional form, Eq. (2.15), the aircraft motions 
were assumed to be slowly varying. The nonlinear indicial response for 
these arbitrary large-amplitude motions was assumed to be the same as the 
response to a motion with fixed past. For flutter motions involving 
small-amplitude, high-frequency oscillations of the motion variables about 
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fixed mean values, the assumption of a constant past history is also 
justified. Here the maximum excursions of the motion variables are small 
and the motion can be considered to be generated by a series of steps 
applied not at the instantaneous values but rather at the mean values of 
the motion variables. The approximate integral form is thus valid for 
flutter motions, but cannot be simplified as was done in the case of slowly 
varying motions. 

When considering large-amplitude motions of higher frequencies, the 

simple assumption that the general nonlinear indicial response is the same 

as the response for fixed past may no longer be adequate. A more precise 

way to account for the dependence of the indicial response on past events, 

for the planar motion discussed above, is to assume that the response to a 

step at time t is dependent not only on the levels cx(t), V(i), and 

q (t), but also on their rates of change at the time of the step, a(r), 

V(t), and q (t). This is illustrated in Fig. 4 for a step in a. In 
B 

motion C, a, V, and vary linearly in time with rates S(t), V(t), and 

Jo 

q (t) over the time interval i-t < £ < t during which events may influ- 
enee the subsequent response. The response to motion C is assumed to 
be a closer representation of the response to the arbitrary motion A 
than is the previously discussed response to motion B, whose past history 
is held fixed. The approximate response is again a function rather than a 
functional and is dependent on the parameters a(i'), a(x), V(t), V(t), 
q-ft), and q (t). Since the motion is uniquely specified by these param- 
eters over the interval of influence prior to the step* the response is 
again dependent on the time variable t-T rather than on t and 1 
separately. In functional notation this is indicated as 
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Ctt A [o(e),V(?),q. (5);t,T] 

S 

« C m ^[a(r)+a(T) (5 -t) ,V(t)+V(t) (C-T)^ q B (T)+4 B (t) (C-x) ;t,x] 

A a 

= C m ^(a(r) ,a(t) ,V(x) ,V( t) , q B (T)4 B (t) ;t-T ) (2.22) 

a 

In addition to more accurately describing large-amplitude, high-frequency 
motions, the expanded approximate indicial response can describe motions 
involving hysteresis effects, where, for example, the response to a step 
imposed at a given level of the motion variables may be different, depend- 
ing on whether the variable was increasing or decreasing prior to the 
application of the step. The expanded approximate integral form corre- 
sponding to Eq. (2,17) is 


c (t> 

m 


I 


= c m (°) + 1 ^ (a(x) ,a(T),V(x) 5 V(x) ^(x) ^(x) ;t-x ) (a)dx 


o a 


B *^B 


dx 


■j: 


+ | 5m v (S(t),a(T) ,V(x) ,V(t) , q B (r) 3 4 b (t) 5 t-x) — (v/V R )dT 

o 


i; 


+ \ ^ (a(x) ,«(x ) ,V(x) ,V( x ) »q. B ( x ) j4 b ( t ) »t-T ) ^ (q. B «,/V R )dx 

o B 


(2.23) 


which is seen to be identical to the simplified form if the indicial 

response functions are found to "be independent of the rates a, V, and 4 . 

B 


2.3 Development for General Body in Free Flight 
One advantage of using the theory of functionals to develop nonlinear 
aerodynamic formulations is the ease with which the analysis may "be 
extended to include additional independent variables. Once the functional 
dependence of the force and moment on the variables is established s the 
succeeding analysis corresponds formally to that presented for the case of 
planar motion. We now use this formalism to develop a nonlinear formulation 
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for a general "body in free flight. Under the same restrictions on the 
altitude and flight speed variations as were imposed on the planar case, 
the force and moment acting on the body are dependent solely on the veloc- 
ity and angular velocity history of the motion. The reactions and the 
motion variables can he expressed in terms of their components resolved 
in the body- fixed axes or, equivalently, by their components resolved in 
the aerodynamic axis system. The resulting nonlinear formulations are 
equivalent, hut lead to different sets of characteristic motions from 
which the total reactions are determined. However, as will he seen, the 
importance of coning motion is evidenced by the fact that it appears as 
one of the characteristic motions of the formulations developed in both 
axis systems. 

2.3.1 Body-Fixed Axes 

The components of the flight velocity vector resolved in the body- 
fixed axes, Ug, v^, w^, are related to a, 6, V through Eqs. (2.1) and 
(2.7)* Thus the expanded dependence of the pitching-moment coefficient 
can be specified as a functional of the form 

C m (t) = E[a( 0 ,£U),VU),p B U),q B U),r B U)] (2.2U) 

The formulation of the nonlinear indicial responses and of the exact and 

approximate integral forms for C^(t) parallels that of the planar case, 

Eqs. (2.lU), (2.15), and (2.17), respectively. When the integral form is 

expanded about the point a = a(t), 3 = £(t) , V = V(t ) , p_ = q^ = r = 0, 

b b b 

and only terms linear in the rates are retained, the resulting nonlinear 
formulation corresponding to Eq. (2.20) is 
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= C m (“;a(t),B(t),V(t))+ vnju 

R 


PT,(t)jl ^ 

C™ («;a(t),&(t).V(t)) 


q. B (t)£ 


r B (t)£ 


V lu a 
R q B 


C m (®;a(t) ,B(t) ,V(t) ) + — — C™ (°°;a(t) ,3(t) ,V(t) ) 

Q.-n V t-. 


R 


'B 


• • 

+ sMi C m .(a(t),g(t),V(t)) + %^C m .(a(t),B(t),V(t)) 

" Ti » T-i O 


R a 


+ cL.(o(t),e(t),v(t)) 

V 

R 


(2.25) 


where, for brevity, the zeros corresponding to p^, q^, have been 
omitted* 

Just as was done in the planar case, the formulation can be further 
simplified for motions where the plunging of the mass center is small and 
where the flight speed makes only small oscillations about a mean speed. 

As before, we rewrite Eq. (2*25) and, guided by Eq. (2*ll), neglect terms 
multiplied by ), yr_+(3-ap_), and V (their contributions vanish 

identically in the case of zero plunging and constant flight speed)* The 
reference velocity V is taken as the mean flight velocity V(t) and is 
omitted from the functional notation of Eq* (2*25) for conciseness* The 
simplified nonlinear formulation is 


Cjt) = ^(“jS.8) + (Cj^ (»;«,£) + yC m .(S,g)) 

B a 


i 3* t% 


^ , , 


Y V wm r 


(Caw. - Y C m ,(« s 6)) + -f- (yCmp 


B 


3 


/V A 

0C, 


nu 


+ oiC m (°°;a,B)) 


£ B 


B 


(2.26) 


Equation (2.26) is identical, to Eq. (22) of Tobak and Schiff, 10 derived 
for the case of constant flight speed. Analogous expressions for the 
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A A A 

other moment coefficients C and C. , and for the force coefficients C v% 

n l s X 5 

C Y , C z are obtained by substituting them wherever appears in 

Eq. (2.26). 

Each term of Eq, (2.26) is associated with a particular motion from 
which it may be evaluated. The first two terms are identified by compar- 
ing them with those previously obtained for the case of planar motion, 

A 

where = 0, 8 = const = 0. The first term is thus the pitching-moment 

coefficient that would be measured in a steady planar motion with a, 8 

(and V) held fixed. The term (C^ +yC m . ) is the planar damping- in-pitch 

q B a 

coefficient that would be measured for small angular oscillations in a 

about fixed a, with 8 held fixed at 8(t) and = 0, Similarly, 

(C m -TC m . ) is the change in the pitching-moment coefficient due to 
r S 3 

damping-in-yaw motion (small angular oscillations in 6 about fixed g, 
with a held fixed, = 0). As was pointed out in Ref. 10, the term 
(C m -lC m . ) and the analogous term in C (t), (C n +YC n .), are the cross- 

b $ n % a 

coupling terms normally excluded in the classical treatment. These terms 

are missed by attempts to generalize to the nonlinear case from linear 

formulations based on the principle of superposition. 

The last term in Eq. (2.26) is identified by comparing it to the 

result that would be obtained for the case of steady coning motion 

(a = const, $ = const, <J> = const), where, as seen from Eq. (2.12), p = yh 

Jd 

A * ^ * 

q_ B = B<j>, r £ = otcfi . When these conditions are substituted in Eq. (2.25), 
the result is 


= 5 («;£,&) + (yC™ (~;a,e) + £c (“;a,g) + aC m (»;a,£)) 

P B q B B 


( 2 . 27 ) 


24 




• » 

<f>, evaluated at <j> = 0, of the pitching-moment coefficient measured in 

steady coning motion, and is designated Cm (°°;a,g}. 

♦ 

In summary, Eq. (2.26) suggests that for the free flight of a general 
body with small plunging and near-constant flight speed, the total moment 
may be compounded of the contributions from four characteristic motions: 
steady angle of attack and sideslip, planar pitch and yaw oscillations at 
constant angles of attack and sideslip, and coning at steady angle of 
attack and sideslip. These motions are illustrated schematically in 
Fig. 5. 

2.3.2 Aerodynamic Axes 

As can be seen in section 2.1, the velocity vector of a general 
motion can be specified in the aerodynamic axes by the scalar variables 
6, V. The angular velocity vector is specified by p , q, r or, 
equivalently, by X, q, r, since p^ is related to x and through 
Eq. (2.1*). The pitching-moment coefficient in the aerodynamic axes, 

C (t), thus may be specified as a functional of the form 

C m (t) = F[6U),*U),vU),xU),<l(s)»r(s)] (2.28) 

Proceeding formally, one finds that the nonlinear formulation in the aero- 
dynamic axes corresponding to Eq. (2.25) is 


C (t) = C (»;6(t),*(t),V(t)) + 
mm V R 

+ C n ^(»;6(t) > 4.(t),V(t)) 

R 

+ ^T^ c m ( <5 (t ) ,ifi(t ) ,V(t ) ) + 
R 6 

+ C m -(fi(t),iMt),V(t)) 

V 

R 


c m. (°°;s(t) ,tfi(t) ,v(t) ) 

X 

+ Cm (-;6(t),^(t),V(t)) 

R 

^%^C m .(6(t),^(t),v(t)) 

R t 

(2.29) 
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As before, Eq. (2.29) may be simplified in the ease of small plunging and 
near-constant flight speed by rearranging terms and neglecting those terms 
multiplied by q-a, r-e\, and V (which vanish identically for zero plunging 
and constant flight speed, e.g. , Eq. (2.10)) to obtain 


• * 

C m (t) = + T ( c m a (“^ 4 t) + yCmJS,*)) + ^ C^U,*) 

q 

IX!, 

+ Y T + 5C m (~;6,1/|)) (2.30) 

Analogous expressions hold for thq, other force and moment coefficients. 

The characteristic motions from which the terms of Eq. (2.30) are evaluated 

differ from those in the body- fixed axes. - Here the term (Cm +YC m ) is the 

^ 5 

planar damping-in-pitch coefficient that would be measured in a nonplung- 
ing motion for small oscillations of the resultant angle of attack a 
about fixed a with tp and V held fixed, X = 0. This term is designated 

c m. ( *5 >V*) • The term is the change in the pitching-moment coefficient 

o ^ 

due to damping-in-roll motion (small oscillations in ip about fixed ip 

with o = const, V = const, X = 0). When the conditions of steady coning 

« • • 

motion (q = 0, r = 6$, A = -y$) are substituted in Eq. (2.29), it can be 

shown that the term (<SC m +yC m ) is the rate of change with coning rate 

* 

4>, evaluated at <J> = 0, of the pitching-moment coefficient that would 

be measured in coning motion, and this term is designated C m (°°;6,<p). 

* 

Thus 


C m .U,>p) = C m (-6, ip) + YCm.(M) (2.31a) 

a ^ 6 

(“;5 ,ip) = YPni. + 6C m {»;S,*) (2.31b) 

X r 


In the aerodynamic axis system, the characteristic motions are: steady 

resultant angle of attack and bank angle 0, pitch and roll oscillations 
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at steady angles of attack and "bank* and coning at steady resultant angle 
of attack with fixed "bank angle. These motions are illustrated 
schematically in Fig. 6. 

2.3.3 Correspondence Between Axis Systems 

The formulation developed in the body-fixed axis system* Eq. (2.25)* 
is related to the one developed in the aerodynamic axes* Eq. (2.29)* 
through Eq. (2.9)* These relations* presented in detail for convenient 
reference in appendix A* can be used to show the following significant 
equivalence: 

“ Y c n.( 5 >’^ + <$C m> (6,^) = <${(6™ (°°;a,l) + yC m ^(a,B)) 

<fi ^ a B a 

+ (C n (°°;a,8) - Y c n x (a 9 8))} 

B 8 

(2.32) 

The term C n< is the rate of change with $£/V, evaluated at = 0, of 

$ 

the side-moment coefficient C n that would be measured in steady coning 
motion, while the term C n> is the change in the side-moment coefficient 
due to damping-in-roll motion. Thus a determination of C n> -YC n> would 
be equivalent to a determination of the three planar damping coefficients. 
The identity is shown schematically in Fig. 7. There are two cases of 
special interest that lead to simplification of Eq. (2.32): when ^ = 0 
(a = 6, 8 = 0), and when ip = tt/ 2 (a = 0, 8 = $)• In the first case, 

Eq. (2.32) becomes 

C n< (»;6,0) - YC n> («,0) = 5{C n (»;a,0) - yC^U.O)} (2.33a) 
l ^ B 6 

while, in the second case, we obtain 

C n .(®;6,ir/2) - Y C n . (6,ir/2) = 6{C («;0,p) + yC^.^O^)} (2.33b) 

4> ^ “ a 
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A determination of (C n _-yC n> ) as the resultant angle of attach "becomes 
small (6 0) for these cases is seen to he equivalent to a determination 

of the classical (linear formulation) damping-in-yav and damping-in-pitch 
coefficients, respectively. 

Another important equivalence between terms of the formulations in 
the two-axis systems is 


Cm. ( = {Cjn (“iajfJj+yCja. (a, 3) }cos 2 0 +{C n (“;a,3 )-yC nx (a ,6) }sin 2 ij; 


B 


'B 


ru 

& 


-[{C n (“ja,g)+yC n . (a^WC™ (®;a,B J-yC^ci ,g) } ]sin ij; cos 1 ); 


A B 


B 


8 


(2.3*0 

which in the special cases of ip = 0 and = ir /2 becomes, respectively. 


C m , (6,0) = {C m (»;a,0) + yC m .(S,0)} (2.35a) 

a a 

( <5 ,tt/2 ) = (C n («;0,8) - yC .(0,8)} (2.35b) 

a r B 3 

Here, in contrast to Eq. (2.33), the determination of as the angle 

a 

of attack becomes small is equivalent to the determination of the linear 
formulation damping-in-pitch and damping- in-yaw coefficients, respectively. 
2.3.*+ Bodies of Revolution 

When the body has axial symmetry about the x_ axis, the value of 

Jj 

♦ is arbitrary and will be chosen equal to tt /2 for convenience. Also, 

in the case of a body of revolution, the damping-in- roll motion used to 

evaluate Cm. and C n- can be replaced by the classical Magnus experiment 

$ 

where the body is placed at angle of attack and spun about the x axis 
(a = const, ^ = const, X = 0). Here the term C n# is the rate of change 
with , evaluated at ^ = 0 , of the side-moment coefficient that 

would be measured in the Magnus experiment. The only mechanisms for 
affecting the moments on the body in this experiment is the small 
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asymmetry of the flow field produced by viscous shear at the "body surface. 
Thus numerous experimental investigations (e,g, , Regan and Horanoff, 14 
Platou and Sternberg, 15 Platou 15 ) and a viscous theory (Sedney 17 ) have 
shown that, when measured in this manner, the term C n> and the corre- 
sponding side force term, Cy. , for a body of revolution are extremely 
small. Additionally, Schiff and Tobak 5 have demonstrated experimentally 
that, at least for a 10° half-angle cone, the term C n> can be neglected 

'I' 

in comparison with C n> , the change in the side moment due to steady con- 

£ 

ing motion. In this case the determination of C n> alone, for small 
resultant angle of attack, is seen from Eq. (2.33b) to be equivalent to 
the determination of the linear damping-in-pitch coefficient. 

The form that the nonlinear formulation takes when applied to the 
nonplunging flight of a body of revolution illustrates an erroneous assump- 
tion that has frequently been made when attempting to generalize to the 
nonlinear case from linear formulations based on the principle of super- 
position. Here the terms C mil C m> , C n> , and the static side-moment term 

^ o 

C n (“;5,<Ji) which appear in Eq. (2.30 ) and in the analogous expression for 
C n (t) can normally be neglected on the basis of symmetry arguments. With 
chosen as ir/2, the resulting formulation can, with the aid of Eqs. 

(2.10c) and (2.11c), be written as 

ft 

C (t) = C m (»;6 ,tt/2) + — C m> (6,tt/2) (2.36a) 

a 

Xo 

C n (t) = V ^n. - yC n> (6 ) + ~y— C n> (6 ,tt/ 2) (2.36b) 

l ^ i 

Although, as discussed above, C n could be neglected in comparison with 

$ 

C in the first term of Eq. (2.36b), we shall retain it in the notation 
n* 

$ 

to avoid confusion. When Eqs. (2.33b) and (2.35b) are substituted in 
Eq. (2.36), we obtain 
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(2.37a) 


C m (t) = + ^ { Cn r (~50.i) - yC n , (0,6)} 


cjt) = % + yC m .(o,B)} 


c n ( 6 ,tt/ 2 ) (2.3Tb) 

* 


When the linear formulation is valid (i.e., as a 0, g -+ 0), the axial 
symmetry of the body dictates that 

{C m (°° ;0 , 0-K) ) + Y C m (0,6^0)} = (C n (»;0,6-K)) - yC n . (0,g-K))} 


This leads to the linear formulation result that 

C n." ^C n . 

c m. * i - %, t») + 


( 2 . 38 ) 


(2.39) 


where C~ («) + yC^. denotes the planar pitch-damping coefficient evalu- 

% a 

ated at ot = 0, 3 = 0. As the resultant angle of attack becomes large, 
relation (2.38) is no longer valid. The damping measured for pitch oscil- 
lations at large sideslip angle (and a = 0) is not necessarily equal to 
the damping measured for yaw oscillations at the same angle of sideslip. 
Thus, for large 6, need not he equal to (C n- - YC n# )/6 and, in fact, 

a } J 

neither need equal the linear formulation value of the pitch-damping coef- 
ficient. That this inequality does actually occur is demonstrated, in Fig. 

8, which compares the values of SC^ measured by Iyengar 4 for a 10° 

o 

half-angle cone at Mach number = 2 to the values of C n ^ measured in 

4 > 

steady coning motion at the same conditions, where C n- << C n# (Kef. 5). 

- A ^ 

Also shown is the theoretical value of 6{C m (») + y6 m .} (Tobak and 

S 

Wehrend ) . At the low values of 6 the three values are in good agree- 
ment. However, as 6 increases the equality breaks down. Consequently, 
extensions to the linear formulation that retain the equality between 0^ 

a 

and the term (C n> - YC n- )/<$, although they allow for the nonlinear 
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"behavior of C with increasing angle of attack, will lead to false 
a 

results* In particular, when (C n< - yC n> )/6 is incorrectly forced to 

4> Ji 

equal C m , the difference between the assigned value and the actual one 
c 

must be absorbed in the remaining term of the side-moment equation, 
{p^fc/V}C n<6 As discussed by Levy and Tobak,^ this may cause methods that 
extract nonlinear aerodynamic coefficients from free-flight data, if based 

on such an aerodynamic formulation, to assign an unrealistically large 

\ 

value to C n , although the value determined in the classical Magnus 

\ ' \ 

experiment 1 would b l e negligibly small* 

\ \ \ 

2.3^5 Potential Application to Airplane Spins 

i ■ 

The ability to predict the pre-stall and post-stall spin behavior of 
high perform l ance aircraft is currently-" hampered by the inadequacies of the 
linear aerodynamic moment system. The striking similarity between coning 
motion and the steady spin of an aircraft suggests that a moment formula- 
tion similar to Eq. (2.25) (or, alternately * Eq. (2.29)5 will properly 
describe the aerodynamic reactions on a spinning airplane. It is known, 
however, that in the establishment of a spin the large asymmetric regions 
of separated flow on the wings of the vehicle cause the aerodynamic reac- 
tions to be highly nonlinear functions of the spin rate, even at low spin 
rates. This contradicts the assumption, used in the development of Eq. 
(2.25), that' the reactions are linear functions of the rates. It is 
anticipated that by expanding the integral form for C^Ct) corresponding 
to Eq, (2.17) about the point p g = Pg(t) (the instantaneous spin rate) 
rather than about p = 0, a formulation corresponding to Eq. (2.25) could 

ID 

be obtained which would describe the nonlinear behavior of the aerodynamic 
reactions with coning rate as well as with angle of attack and sideslip. 

In such a formulation p^Jt) would, of necessity, be retained in the 

' .D 

notation . 
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3. NUMERICAL FLOW-FIELD SOLUTION 


The previous chapter has shown that the nonlinear aerodynamic force 
and moment acting on a body performing large-amplitude nonplanar motions 
can be compounded of the contributions from four characteristic motions: 
steady angle of attack* pitch oscillations* either roll or yaw oscilla- 
tions* and coning motion. It would be desirable to be able to apply 
methods of computational fluid dynamics to compute the flow fields about 
bodies performing these characteristic motions and so obtain their con- 
tributions to the aerodynamic reactions. Extensive study of the steady 
angle-of-attack case has led to the development of many numerical finite- 
difference methods for computing the steady inviscid flow field about two- 
and three-dimensional bodies for subsonic through hypersonic Mach numbers. 
The computation of the nonsteady flow fields generated by the oscillatory 
motions is more difficult. Here the flow variables are functions of time 
as well as of position, and thus the solution requires increased computa- 
tion and larger computer data storage capacity than is necessary for the 
steady flow case. Coning motion* shown to have special significance in 
the nonlinear moment formulation, generates a flow field more amenable to 
numerical solution than do those of the oscillatory cases, since to an 
observer fixed on the coning body the surrounding flow is steady. Hence* 
techniques developed for the solution of steady flow fields can be applied 
to coning motion. 

The study of steady supersonic flow has led to the development of a 
class of accurate numerical finite-difference methods termed inarching 
methods. By taking advantage of the fact that the equations governing a 
supersonic flow are hyperbolic in the streamwise direction, the use of 
such methods advances an initial solution, specified at one transverse 
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plane, in the streamvise direction to obtain the entire flow field. These 
methods enable the solution of the full nonlinear gasdynamic equations, 
rather than the simplified equations obtained by the introduction of a 
velocity potential. The solutions obtained are therefore valid for flows 
at hypersonic Mach numbers where vorticity is generated by strong curved 
shock waves, as veil as for flow at lower supersonic Mach numbers. One 
such method, the noncentered second-order scheme introduced by MacCormack 1 1 
and developed by Kutler and Lomax 12 as a shock- capturing technique, has 
been shown to be both accurate and versatile. Results of computations of 
the complex steady hypersonic flow field surrounding a proposed space 
shuttle orbiter, obtained using this technique, show excellent agreement 
with experiment and with results obtained from a three-dimensional method 
of characteristics (Rakich and Kutler 19 ). In a shock-capturing technique, 
the equations are expressed in conservation-law form and the finite- 
difference scheme is applied uniformly at all points of a computational 
mesh which extends into the undisturbed free stream ahead of the bow shock 
wave* The jump in the flow variables across the shock is spread over 
several points of the mesh. In a sharp-shock technique, the bow shock is 
treated as a discontinuity and the Rankine-Hugoniot relations are used to 
determine the flow conditions immediately behind the shock. One advantage 
of the shock-capturing technique is its ability to determine the position 
and strength of the bow shock without special computer coding. A second 
and more important advantage is its ability to determine the position and 
strength of embedded shock waves , such as the crossflow shocks that are 
seen on the leeward side of a body at large angles of attack, if they 
occur within the flow field. Because of its accuracy and simplicity, the 
shock-capturing technique was extended to the case of a body in coning 
motion in a supersonic stream. 
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3-1 Method of Solution 


The nonlinear Eulerian gasdynamic equations are solved numerically 
to determine the totally supersonic inviscid flow field about a pointed 
body in coning motion- A body-fixed cylindrical coordinate system, desig- 
nated the computational system, is introduced. The origin of these coordi- 
nates lies at the center of gravity of the body, with s, the axial coordi- 
nate , aligned with the negative axis. The radial coordinate t lies 

in planes normal to the axis (i.e., in crossflow planes) as illustrated 

in Fig- l(a). The circumferential angle, 9, is seen to be measured from 
the resultant angle-of-attack plane. In coning motion at fixed coning rate 
5>, the flow is time-invariant with respect to an observer fixed in the com- 
putational coordinate system. Since the flow field is everywhere super- 
sonic, the gas dynamic equations are hyperbolic in the axial direction. A 
computational mesh is established between the body surface and the free 
stream ahead of the bow shock wave in planes normal to the s axis. With 
the flow field specified at an initial data plane, s = initial* ^ ac ^ ormac ^ !s 
method is used to inarch the solution in the s direction over the length 
of the body. At the outer edge of the mesh, the flow variables are assigned 
free-stream values, while at the inner edge, the flow is kept tangent to 
the body. With the complete flow field thus determined, the forces and 
moments are obtained from a subsequent integration of the surface pressure 
distribution. 

Computations have been carried out for conical bodies of circular and 
elliptical cross section. An approximate initial solution was generated 
at the initial data plane by assuming the flow upstream to be that about a 
cone at angle of attack and yaw, with a uniform sidevash velocity. Note 
that although the body geometry is conical, the free-stream sidewash 
generated by coning motion is a function of axial position along the body, 



and the resulting flow field is not conical. Details of the method and a 
discussion of the approximate starting solution are given below. 

The computations were carried out on an IBM 360/67 computer linked to 
a cathode-ray tube graphics device. The graphics unit, which possesses 
man-machine interaction capability, was used to study the flow field as it 
developed, and to control any numerical instabilities that evolved. 
Approximately 50 minutes of computer time were required per case (one 
angle of attack at one coning rate). 


3.2 Gasdynamic Equations 

The equations governing the -unsteady inviscid flow of a non-heat- 
conducting perfect gas around a body performing an arbitrary motion, writ- 
ten with respect to a body-fixed coordinate system whose origin is at the 
center of gravity of the body, can be expressed as 

( mass ) 3 ! (p) + div{p?) = 0 (3.1a) 

(momentum) 

(pv) + grad p + P V • grad V + V div(pV) + p j^2QxV + (V cg ) 

+ (q)xV + 5x(^xr)J = 0 (3.1b) 

(energy) 


(£) + div[(p+&)v] + (v cg ) + -g- (Q)xV + nx(Hx r )J . p v = 0 

( 3 . 1 c) 

(state) p = (y-l)pe ( 3 . Id) 

where 


r position vector in the body coordinate system 

V(r,t) velocity vector (having components u, v, w) of the fluid at the 

point r, measured relative to the body-fixed coordinates 
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v Q velocity and angular velocity, respectively, of the 

eg* 

body-fixed coordinates measured with respect to 


an inertial system 

g. = p[e+(l/2) |v| 2 ] total energy per unit volume of fluid 


The terms 2pQxy and pftx(fix?) in Eqs. (3.1b) and (3.1c) are Coriolis and 
centrifugal force terms that appear because the body- fixed coordinate sys- 
tem is noninertial. 

In the case of steady coning motion, V = const, = const, and the 

eg 

flow field is time-invariant in the body- fixed coordinates. Under these 


conditions, Eqs. (3.1) become 

(mass) div(pV) = 0 (3.2a) 

(moment cun) 

grad p + pV • grad V + V div(pV) + p[2fixV + ftx(ftxr)] = 0 (3.2b) 


(energy) 


div[(p+S)Y] + [2ftxV + ftx(ftxr)] ♦ pV = 0 


(3.2c) 


Equation (3.2c) can be simplified upon recognizing that (fixV) • V = 

. (7xV) = 0. Also it can be shown that with ft = const, 
curl[ftx(ftxr) ] = 0, and thus ft*(ftxr) can he expressed as the gradient of 
a potential, i.e., ftx(ft x ?) = grad(-$/2). Substituting in Eq. (3.2c), and 
using Eqs. (3.2a) and (3. Id), we obtain 


This states that 


pV * grad^ye -•|‘£+-|'lv| 2 ^=0 
(ve - I ♦ + f M 2 ) is constant along streamlines 


(3.3) 


of the 

flow field. Since the motion under consideration is that of a tody through 
a uniform atmosphere at rest with respect to the inertial system, the con- 
stant is the same for all streamlines, and Eq. (3.3) "becomes 
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1 1 

ye - j V { 2 = const = ye^ = h Q ( 3 . ha) 

or 

p = p[ 2 h Q + <f - (u 2 +v 2 +w 2 )] ( 3 .hb) 

The number of dependent variables is thus reduced from five (E, p, V) to 
four (p, V) , with p related to p and V through Eq. (3.hb). 

When the coordinate system is specialized to the body-fixed computa- 
tional system, the gasdynamic equations ( 3 . 2 a) and ( 3 . 2 b) can be written 
in conservation-law form as 

E 1 + F’ + G’ + H' = 0 (3.5) 

S T O 

where the subscripts denote differentiation and 


E’ = 

pu 

p+pu 2 

jpt = 

pv 

puv 

, G T = — 

pw 

pw 


puv 

7 

p+pv 7 

■ T 

pw 


puw 


pw 

1 

p+pV 2 


pv 

puv* + pr[2(aj2V-ai3v)4'ai 1 ai2T-s (t02 + ^3 ) 1 
p(v 2 -w 2 ) + pt[ 2 ((j 03 U-^ 1 v)+^ 1 o) 2 S-t(( j ji+o) 3 ) ] 

2pw + pt[2(u)iy-lo 2 u ) +w 3 ( u) 2 r+to l s ) ] 

The components of the angular velocity vector of magnitude resolved 

* 

in the s,x ,0 directions are, respectively, u>i = -<$> cos a, 

* • 

a >2 = sin a cos 0, and 0)3 = -<p sin a sin 8 . The energy equation is given 
by Eq. (3.hb), while the centrifugal force potential, obtained from the 
integration of grad(-$/ 2 ) - nx(Hxr), can be expressed as 

$ = $ 2 [(s sin a + x^cos a cos 0) 2 + (t sin 6 ) 2 ] 

For supersonic flow, Eq. (3.5) is hyperbolic with respect to s, and 
thus can be integrated in the s direction. Specifically, u, the compo- 
nent of the local flow velocity in the s direction, must be greater than 
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the local speed of sound at all points of the flow field. In regions 
where this condition does not hold, such as in the subsonic nose region of 
a blunted body, the marching method fails. In such a case the nonsteady 
form of the gasdynaraic equations would have to be integrated with respect 
to time, and the steady solution would be obtained as the steady limit of 
an unsteady flow. 

It is generally advantageous, when flow problems are solved with the 
use of finite-difference methods, to have the physical boundaries of the 
flow field lie along coordinate surfaces of the computational coordinate 
system. Under these conditions the application of the surface boundary 
conditions is greatly simplified. Further, to improve the accuracy of the 
solution, it is desirable to use a dense spacing of computational grid 
points in those regions of the flow field where large gradients of the flow 
variables are known or suspected to exist. Thus, in the present case, that 
of flow over cones of circular and elliptical cross section, the annular 
region of interest about the body (in the t,0 plane) is transformed into 
a rectangular region (Fig. 9). A radial independent coordinate y is 
chosen to map the region between the body surface and an outer boundary 
into the region 0 < y < 1. The outer boundary is chosen to lie in the 
undisturbed stream outside the bow shock wave. It is known that when a 
cone of elliptical cross section is placed at incidence in a wing-like 
attitude (i.e. , with the semimajor axis of the ellipse normal to the flow 
direction) rapid variations of the flow variables occur in the vicinity 
of the semiraajor axis. A circumferential independent variable y, depen- 
dent on the body cross section, is therefore chosen to cluster the physi- 
cal circumferential planes more closely in these regions. The 
transformations are: 
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S S 9 S initial “ S S final 

(3.6a) 


y = (x-T b )/(T ob -T b ) ; 0 * v * 1 

(3.6b) 

where 

n = tan -1 ^| tan ej ; 0<n<2ir 

( 3 . 6 c) 

T b = T h {e ’ s) 

value of t at the "body surface 


T oh = T ob (0 ‘ s) 

value of t at the outer boundary 


A 

sine of the angle subtended by the semimajor 
of the elliptical cone 

axis, a(s) 

B 

sine of the angle subtended by the semiminor 
of the elliptical cone 

axis, b(s) 


Applying these transformations to Eq. (3.5) gives 

E + F + G + H = 0 (3.7) 

s y n 

where 
E = E' 


F « 7-I7- j[(v-l)T^ s -yT 0b ]E’ + F- + [(y-l)r b -yr ob ]G',' 
ob b ( 9 


cos n | 

' A s 


, A 2 cos 2 n + B 2 sin 2 r) 


b) s 

+ AB G 

|cos 2p 1 

( B s 

A aV 

T 0bo“ T b s l r A 2 -B 2 

T ob- T b J E ' + l Sln 2n » 

u 

--) + 


T ohQ- T he 1 
T oh“ T h -J 


G* 


Here denotes differentiation of B with respect to s, x ob0 denotes 

differentiation of t ^ with respect to 6, etc. 


3.3 Differencing Scheme 

The rectangular region of interest in the y,n plane is divided into 
an equispaeed rectangular grid having 20 intervals Ay in the y direction 
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and 36 intervals An in the n direction. Thus Ay = 0.05, and 
An = tt/18. With the primitive flow variables (pressure, density, and 
velocity) known at all mesh points of an initial data plane at s * n 4s, 
MacCormack 1 s two-step, predictor-corrector difference scheme is used to 
advance the solution to s = (n+l)As. As applied to Eg.. (3.7), this method 
is 


s n+ i 


= E' 


,n 




— /p n - F n ^ - — (& n 

Ay\ j+l,k j»k / An \ j,k+l 



AsH 1 ? 

J » k 


(3.8a) 



1 

2 



+ E n+1 


j > k 


As /^n+l «n+l \ As_/~n+i _ ~n+l \ 

AU V j,k ' Vi, k ; ~ An ,k J, k-i/ 


. „n+l 
AsH. , 
J , k 


(3.8b) 

where E 1 } = E(n As , j Ay ,k An), the predicted value of E at s = n As , 

J 5^ 

E 1 ) = E(n As,;) Ay,k An), the corrected value of E at s = n As , 

j 

F 1 ? = F(E^ ) , F? = F(E^ ), etc. To proceed, the conservation-law 

j ,k j ,k j ,*£ j ,k 

variables E 1 ? . , F 1 ? , etc., are formed from the primitive variables at 

<3 »k J ,k 

all points of the mesh since the body and outer boundary geometry are 

-ft* 1 

known. Predicted values E. , of the conservation-law variables at 

J ,k 

s = (n+l)As are obtained with the use of Eg. (3.8a), and these are decoded 

to obtain the predicted values of the primitive flow variables. Boundary 

conditions are then applied to ensure tangency of the flow at the body, 

and the proper flow conditions in the free stream. The conservation-law 

variables , G-f. 1 , and are formed from the predicted primitive 

variables and the known body and outer boundary geometry at the new axial 

position. Corrected values E 3 ? 4 ? are obtained using Eq. (3.8b), and are 

J 

decoded to obtain the solution in terms of primitive variables at 
s = (n+l)As. The boundary conditions are again applied, thus completing 
the computational cycle. This cycle is repeated to advance the solution 
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from s = (n+l)As to s = (n+2)As 4 and so on over the entire length of 
the body. 

The primitive flow variables are obtained from the new (predicted or 
corrected) values of E and the energy equation (3.tb) in the following 
manner. Let 


pu 


e l 

p+pu 2 

= 

e 2 

puv 


S3 

puw 


S4 


(3.9a) 
(3. 9b) 
(3.9c) 
(3.9d) 


With pressure eliminated from Eq. (3.9) by the use of Eq. (3.^b) 4 the 
simultaneous solution of Eq. (3.9) gives 


v = e 3 /e! 
w = eit/ei 

2 

u = JL (7_x_ izi 

Y+l |_V +1 e V y+l 

p = ei/u 


(3.10a) 

(3.10b) 

(2h Q + $ _ v 2 - w 2 )j (3.10c) 

(3.10d) 


The pressure p is determined from the energy equation (3.H). The posi- 
tive sign appearing before the radical in Eq. (3.10e) is used because the 
axial component of the local flow velocity is supersonic throughout the 
flow field. That the flow was supersonic in the cases considered was 
assured since the maximum effective deflection angle, defined as the sum 
of the resultant angle of attack and the cone half-angle, was always less 
than the shock detachment angle for axisymmetric supersonic cone flow at 
the same Mach number. 

It is desirable, both for reasons of numerical accuracy ,and for com- 
putational efficiency, to take as large a step size As as possible. 
However, explicit numerical schemes such as this one become unstable if 
too large a step is employed. To determine the maximum step size that may 



"be used at a particular s station, the nonlinear equation (3.7) is 
expressed in a locally linear form as 


E + ME + HE + H = 0 (3.11) 

sun 

where M and N are the Jacobian matrices 3P/3E and 9G-/9E, respectively. 
The maximum allowable step sizes determined from applying linear one- 
dimensional stability analyses successively in the y,s and n,s planes 
(e.g., Richtrayer and Morton 20 ) are 

(3.12a) 


As 1 

* I a M L 


As 


M 1 max 
1 


An ' JaJ 


(3.12b) 


N 'max 


where and a ^ are the eigenvalues of M and N, The eigenvalues are 
evaluated for all points in the flow field at the current s location, 
and the most restrictive condition on As is selected* While the use of 
uncoupled one -dimensional analyses is not exact, the results are useful 
provided the step size thus determined is further reduced by approximately 
10 percent* Further, if the local Mach number tends toward unity, the 
step size determined in this manner tends toward zero, thus providing a 
warning that the assumptions of the analysis may no longer be valid. 


3.^ Boundary Conditions 

The outer edge, y = 1, of the computational grid is chosen to lie in 
the undisturbed free stream ahead of the bow shock. Here the pressure and 
density are the free-stream values, while the velocity components, deter- 
mined from geometrical considerations , are 

u = V cos a - 4>t , sin a sin 0 
ob 

v = -V sin a cos G + 4>s sin a sin 6 
w = V sin a sin 0 + $(s sin a cos 0 + cos a) 


h2 


(3.13a) 

(3.13b) 

(3.13c) 



At the sides (n = 0 and n = 2 tt) of the computational grid, a periodic con- 
tinuation principle is applied, i.e., p(y,r) = 0) - p(y,n = 2tt), 
p(y,0 = 0) - p(y ,n =s 2 tt), etc. The tangency boundary condition at the 
body, y = 0, is due to a scheme of Abbett 21 - and is briefly summarized here. 
The flow variables are known at the new axial position after the predictor 
step. In general, the local flow velocity at the body is not tangent to 

the body. A local two-dimensional Pr an dt 1 -Mayer expansion or an isentropic 

* 

compression is used as needed to turn the flow into the local tangent plane. 
This satisfies the tangency condition and determines a corrected value of 
the surface pressure. A corrected value of the surface density is then 
determined since it is assumed that flow along surface streamlines is 
isentropic. The magnitude of the corrected flow velocity is then deter- 
mined from the energy equation, Eq. (3.4), while its direction is that 
obtained from satisfying the tangency condition. This tangency condition 
scheme has been shown to give good results. The surface flow conditions 
for a circular cone at angle of attack computed by Kutler et at. 22 using 
this scheme were found to be in good agreement with the numerical results 
of Babenko et al. 22 and with the method of characteristics results of 
Rakich. 24 

Crossflow shocks were observed to exist in the flow field for the high 
angle-of-attack computations and to extend to the body surface. In such 
cases, the assumption that the flow is isentropic along surface streamlines 
is no longer valid. In principle, the value of surface entropy (used in 
the tangency scheme to compute the corrected surface density) for points 
behind the embedded shocks should be different from those points unaffected 
by the shocks . However, if the embedded shocks are relatively weak this 
refinement to the tangency condition can be justifiably neglected. This 


43 



was found to "be the case in the present computations, and a single value 
of surface entropy was used for all points on the surface. 


3.5 Initial Solution 

The flow field surrounding even a conical body is not conical when 
the body is in coning motion since the sidewash velocity of the undisturbed 
stream, measured by an observer on the body, is a function of axial position. 
This is illustrated in Fig. 10(a) for the conical tip of a body in coning 
motion. The velocity component normal to the angle-of-attack plane varies 


linearly from <Jis. . sin a at s = s , . to <f>s. ... , sin a at s = s. ... . , 
J tip tip initial initial’ 

changes direction at s = 0, and increases linearly with s behind the 
mass center. The exact computation of the flow field upstream of the ini- 
tial data plane at s = would require the use of a three- 

dimensional, time-dependent technique with the steady solution approached 
as the limit of an unsteady flow. To avoid this complication, an approxi- 
mate initial solution, valid for bodies with conical tips and capable of 
being generated by the inarching method, is employed. The flow field 
upstream of the initial data plane for the conical tip in coning motion 
(<t> = const) is assumed to be that of the body in steady planar motion at 
the same values of angle of attack and bank, with a uniform imposed side- 


wash velocity of ^ s j_ n pti a i s ^ n a ac_t ing normal to the angle-of-attack 
plane. This is illustrated in Fig. 10(b). Under these conditions the 
flow field around the conical tip is conical. 

The marching method is used to generate the approximate conical solu- 
tion in a "distance asymptotic" technique. The outer boundary of the flow 
field is chosen as a cone whose apex is coincident with that of the body. 
Corresponding points in successive s = const planes thus lie along rays 
of the flow field. The flow variables are set to convenient values, such 



as free-stream values , and the gasdynamic equations are integrated down- 
stream. As the integration proceeds, the arbitrarily chosen initial 
values of the flow variables have less and less effect on the current 
values. The current solution is increasingly affected only by the 
boundary conditions at the body surface and at the outer boundary. When, 
at corresponding points, no change in the flow variables is detected with 
further integration downstream, the flow variables are constant along rays 
and the solution has been determined. The conical field is then scaled 


to place it at s = s. ... ... 

initial 

’When the surface tangency condition is applied during the determina- 
tion of the conical starting solution, the value of the surface entropy 
is allowed to vary. At each step the position of the crossflow stagnation 
streamline is determined numerically, and the surface entropy is chosen 
as that of stagnation streamline one mesh interval above the body. As 
the conical solution is approached, the surface entropy, too, remains 
constant with further integration downstream. The numerical method dif- 
ferences across the vortical layer known to exist in the flow fields about 
circular and noncircular cones at angle of attack (cf. Ferri 25 ). Although 
differencing across these regions of high gradients causes some small 
spatial oscillations of the flow variables close to the body, it has no 
effect on the stability of the method. 

The errors associated with the use of the approximate initial solu- 


tion should he small, particularly when initial chosen close to . 
s t’p* an< ^ siloulci diminish as the integration (for the body in coning 
motion) proceeds downstream from the initial data plane. Further, these 


errors should have only negligible effect on the forces and moments since 


the surface area of the body is small near its apex. The magnitude of 
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these errors was assessed "by* carrying out a series of computations at one 

coning rate, Mach number, and angle of attack, in which the position of 

the initial data plane was varied . No changes in the forces and moments 

were observed as s* ... _ was varied for values of 

initial 

tip 

thus started 5 percent downstream from the nose for the remainder of the 
cases considered. 


) less than 0.05* The solutions were 


^initial s tip^ s final S 


3*6 Results 

Flow fields were computed for a 10° half-angle circular cone in con- 
ing motion at M - 2 and M - 10, and for an elliptic cone in coning motion 
at M = 2. These cases were selected both to demonstrate various capabil- 
ities of the numerical method and to enable a comparison of the numerical 
results with experimental measurements (and indirectly with the results of 
other analytical and numerical methods). 

3*6,1 10° Circular Cone at M “2; Nonlinear vs. Viscous Effects 

Computations were carried out for a 10° half-angle cone in coning 

motion for angles of attack ranging from 0° to 25° (0 < S < 0.^2) and for 

values of the coning-rate parameter (J>£/V ranging from 0 to 0.15. The 

center of mass (s = 0 ) was located at 6l percent of the body length from 

the nose; thus s, . = -0„6l, s. .. . - = -0.56, and , = 0.39* These 

tip initial ’ final 

conditions duplicate those of the experiment described in Ref. 5* Briefly, 
in that experiment a six-component balance mounted on a rotating sting was 
used to measure the forces and moments acting on a cone in coning motion. 

In addition, the viscosity-induced vortex patterns present on the leeward 
side of the body at high angles of attack were investigated using the vapor 
screen flow visualization technique (Allen and Perkins 26 ). 
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3 .6. 1.1 Flow-Field Results 

3.6.1.1a Surface pressure - Circumferential plots of the surface 
pressure for the conical initial solutions generated at s = = 

-0,5^ are presented in Fig. 11. These plots show the surface pressure , 
normalized hy the free-stream pressure, as a function of the circumferen- 
tial angle 0 for each of the angles of attack a computed. These pres- 
sure plots are not symmetrical about 9 = 0° due to the presence of the 
imposed sidewash velocity, but rather are symmetrical about the direction 
of the resultant crossflow velocity vector, i.e., tan ® syrn “ -is initial^* 
Results for two values of the coning rate $ are presented in the figure. 
A value of <j>/V = 0.12 was used for a < 12.5° , and for a > 15°, i/V , vas 

chosen as 0.15. giving values of 0 of 3.8U° and U.80°, respectively. 

? sym 

At angles of attack less than the cone half-angle, the pressure drops 

smoothly from a maximum at the windward ray (e = 0 ) to a minimum at the 

sym 

leeward ray (0=9 + l80° ) . As the angle of attack is increased "beyond 

sym 

the cone half-angle, the surface pressure is observed to drop to a minimum 
approximately 115° from the windward ray, and then to increase gradually 
toward the leeward ray. At angles of attack of 20° and 25°, the surface 
pressure drops to a minimum located 125° and 135°, respectively, from the 
windward ray. The pressure then undergoes an abrupt jump followed by a 
slight increase as 0 increases toward the leeward ray, thus indicating 
the presence of a pair of crossflow shocks extending from the body surface 
into the flow field. 

3^6. 1.1b Crossflow shocks - The strength and location of the cross- 
flow shocks can also be determined from the circumferential plots of the 
surface crossflow Mach number of the conical initial solutions which are 
presented in Fig. 12. The surface crossflow Mach number, M , defined as 



the surface crossflow velocity, , normalized by the local speed of sound, 
is shown as a function of 0 for angles of attach of 15°, 20° , and 25°. 
These plots are symmetric about the windward ray, whose location, indi- 
cated by M c - 0, is again seen to be displaced to “ ^*8° due to the 

imposed sidewash velocity* For angles of attack up to 15° , the crossflow 
velocity is everywhere subsonic* At a = 15° the surface Mach number 
reaches a maximum of 0.88 located 115° from the windward ray, which corre- 
sponds to the location of the minimum surface pressure, and decreases 
smoothly to zero at the leeward ray (0 = l8'5°). At a = 20° and 25°, the 
crossflow velocity is observed to become supersonic approximately 85° and 
75°, respectively, from the windward ray, and reach a maximum Mach number 
followed by an abrupt drop, indicating the crossflow shock, at values of 
6 which correspond to the locations of the jump in the surface pressure* 

In a conical flow field, the crossflow shocks lie along rays from the apex 
and thus the component of the Mach number normal to these shocks at the 
body surface is just the surface crossflow Mach number, M^. The best 
estimates of the upstream (supersonic) surface Mach numbers, M , and 
shock locations, obtained from the numerical results for o = 20° and 25° » 
are indicated in Fig. 12. The corresponding downstream Mach numbers, M , 
which together with the pressure ratio across the shock are determined 
from M using the normal shock relations (e.g. , Liepmann and Roshko 27 ) 

C\ 

are also shown. The pressure behind the crossflow shocks, determined from 

numerically obtained values of the pressure ahead of the shocks and the 

pressure ratio corresponding to M , is shown in Fig. 11 and is in good 

c 1 

agreement with the numerical results. 

3.6.1.1c Flow- field contour maps - Qualitative characteristics of the 


conical flow field, including the strength and location of the bow shock 
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wave and the changes in the flow field due to increasing angle of attack, 
can be determined from the flow- field cross-section contour maps presented 
in Fig* 13. These maps show contours of constant, equi spaced values of 
the pressure, normalized by the free-stream pressure, in the plane 

S ~ S initial for each of the of attack computed* Also shown in 

Fig. 13 are the body cross section and the selected outer boundary of the 

computational region. The contours are generated hy a contour mapping 

program directly from the numerical data. The contour mapping routine 

generates fewer points along contour segments running from the lower left 

to the upper right of the p,n computational mesh than for other segments* 

This relatively large spacing between contour points causes the more Jagged 

appearance of. the contours in the lower right section of the maps. As 

before, the flow field is observed to be symmetric about 8 due to the 

sym 

imposed si dewash velocity, and the location of the symmetry plane is shown 
in the figure. 

At low angles of attack, the bow shock location, indicated by closely 
spaced contour lines within the computational region, is well defined both 
on the windward and leeward sides of the body. As the angle of attack is 
increased beyond the cone half-angle, the windward side shock strength 
increases. On the leeward side, the bow shock becomes weak, tending toward 
a Mach wave, gnd its location is not easily identified. In these cases the 
shock location is indicated by a contour of p/p^ = 1*005. The small closed 
contours seen in the free stream outside the bow shock at angles of attack 
of 20° and 25° are indicative of overshoots in the numerical data to values 
of p/p ^ greater than 1,10, caused hy third-order dispersive errors in the 
numerical scheme. This is also the cause of the jaggedness of the contours 
in the shock layer close to the shock. Experience with higher— order 



differencing schemes that minimize these effects (Kutler et at. 22 ) ^ as 
well as with sharp shock techniques, indicates that the oscillations have 
negligible effect on the solution in regions away from the immediate 
vicinity of the shock. Thus, for example, the surface pressure distribu- 
tion is not adversely affected by the dispersive effects. 

At angles of attack of 20° and 25°, shocks are seen to extend from 
the body surface into the flow. The crossflow shocks are situated at the 
rear of a region of supersonic crossflow located adjacent to the body and 
surrounded by a larger subsonic crossflow region. This is similar to the 
mixed flow pattern that would be observed on a wing section in a transonic 
flow (M m < l). As in the transonic case, the shock is strongest at the 
body, and its strength decreases with distance from the body surface. The 
crossflow shock terminates in the flow field where it joins the sonic line 
delineating the front of the supersonic region. The locations of the cross- 
flow shocks and of the sonic lines, determined from the numerical data, are 
indicated on the contour maps. At a - 20° the shocks extend about 1/2 
body radius into the flow, while at a = 25 ° they extend about 1 body 
radius and, as discussed previously, are stronger. 

3. 6*1. Id Surface pressure contours - Starting with the conical solu- 


tions at s - s initial 9 £ as ^y nam i c equations were integrated over the 


length of the body. The resulting three-dimensional flow fields are 


observed to vary only slightly and gradually with distance along the body. 


To best demonstrate this variation, contours of constant values of the sur- 


face pressure, normalized by the free-stream pressure, are presented in 
Fig. lh for each angle of attack computed. The contours are plotted on 
the developed body surface, which is unrolled at the leeward (0 = l 80 °) 
ray. For zero coning rate (4> = 0), the flow field surrounding the body 
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would be conical, and the flow quantities, including the surface pressure, 
would be constant along rays* In Fig* lh 9 the surface pressure contours 
would merely be straight lines passing through the apex. Thus the devia- 
tion of the pressure contours from such straight lines is a measure of the 
effect of coning on the flow field. For the small values of the coning- 
rate parameter J>2,/V used in the cases computed, the contours are only 
slightly curved, and no large differences are observed between the low 
angle-of-attack and the high angle-of-attack cases. The contours that 
originate part way along the body indicate a gradual change in the pres- 
sure- with distance on the leeward side of the body* Although the curvature 
of the contours is small, the cumulative effect of the curvature is sig- 
nificant* To obtain a quantitative measure of the effects of coning, the 
surface pressure distributions were integrated to obtain force- and 
moment-coefficient data, 

3.6.1. 2 Force and Moment Coefficients 

3.6.1.2a Normal-force and pitching-moment coefficients - The normal- 

force and pitching-moment coefficients are presented in Figs. 15(a) and 

15(b), respectively. These coefficients were obtained from computations 

at fixed angles of attack for various coning rates. At each angle of 

attack, no change in the normal- force or the pitching-moment coefficient 

was observed as the coning-rate parameter was varied* This result, which 

agrees with the experimental observations, indicates that C 7 and 

* 

( do ;^ 9 ^)» the contributions to the normal-force and pitching-moment 

♦ 

coefficients due to coning motion, are negligible. 

Also shown in Figs. 15(a) and 15(b) are the experimental measurements 
(Schiff and Tobak 5 ) of the normal-force and pitching-moment coefficients, 
respectively. The theoretical values (Tobak and Vehrend 18 ) of the linear 
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normal- force and pitching-moment coefficients, Cg (“) and (°°) , are also 

a a 

presented. At low values of < 5 , the computational results, the experimen- 
tal measurements, and the linear theory values agree veil. As 6 increases 
the computational results and the experimental measurements become nonlinear 
functions of 6, and both are greater than the linear theory values. How- 
ever, the computational results and the experimental measurements remain 
in good agreement. 

3 . 6 . 1 , 2 b Side-fovce cmd side-moment coefficients - At fixed angle of 

attack, the computational results indicate that the side-force and side- 

moment coefficients C„ and C are linear functions of the coning-rate 

Y n 

parameter over the range of <j> investigated. This result, which agrees 
with the experimental findings , demonstrates the adequacy of the moment 
formulation, Eq. ( 2 . 30 ), in which only terms linear in the rates are 
retained. 

Since and are linear functions of normalizing by the 

coning-rate parameter yields C Y> (“;6,t|i) and (°°; 6 t ip ) . Computational 

results for Cy. and C n< are presented in Figs. l6(a) and 16(b), respec- 
<P 4 > 

tively, together with the experimentally measured values. Also shown in 

Figs. 16 (a) and l6(b) are the theoretical values (Tobak and Wehrend 18 ) of 

the equivalent linear planar pitch damping derivatives, <$[6.7 (“O+yCg.] 

q B a 

and 5 [C m (“O+yCm. ] » respectively. These equivalences follow from Eq. 

a 

(2.39) since in an inviscid computation for a body of revolution, where 

there is no viscous shear at the body surface, C n is identically zero. 

$ 

At low values of 6, Cy < and C ni are linear functions of d and there is 

£ l 

excellent agreement among the computational results, the experimental 
measurements , and the linear theory values. As 6 increases, the computa- 
tional results and the experimental measurements become nonlinear functions 


52 



of S and depart significantly from the linear theory results. In con- 
trast, the computational results agree veil with the measurements for 
values of 6 up to 0,34, even as Cy peaks and then tends toward nega- 
tive values. As 6 increases beyond 0.35, the computational results 
diverge from those obtained experimentally. 

3.6.1.2c Vapor screen results - Previous attempts to explain the 
experimentally measured nonlinear behavior of the coefficients with 
increasing angle of attack (Schiff and Tobak, 5 Kuhn et at . 28 ) postulated 
that their deviation from the linear theory values was due to the forma- 
tion of a viscosity-induced leeward-side vortex pattern. In contrast, 
the present inviscid computational results , which neglect such a vortex 
pattern, follow the initial nonlinear behavior of the coefficients. This 
suggests that the initial nonlinearities are caused primarily by inviscid 
effects. An examination of the vapor screen studies of the experiment of 
Ref. 5 tends to support the latter hypothesis. Photographs of the leeward- 
side vortex patterns as they appeared in the vapor screen are presented in 
Fig. IT. The vortex cores are visible as dark spots on a light plane 
normal to the velocity vector. Distinct vortices are not visible in the 
vapor screen for values of 6 less than 0.34 (o < 15°). At 6 = 0.34 
(a = 20°), vortices are visible but are small and lie close to the body 
near the leeward ray. Under these conditions the inviscid flow field 
should be a good approximation of the actual one , and the forces and 
moments derived from the computation should, and do, agree well with the 
experimentally measured values. With further increases in the angle of 
attack, the vortices gradually increase in size, and presumably in 
strength, and lie farther from the body. Here the inviscid flow field is 
a poorer approximation of the actual one, and the computationally derived 
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forces and moments "begin to diverge from those measured experimentally. 
Nonetheless, where the viscosity-induced effects are small, the inviscid 
computations agree well with experiment and demonstrate significantly non- 
linear behavior of the forces and moments. 

3.6.2 10° Circular Cone at M = 10 

A computation was carried out for a 10° half-angle cone in coning 
motion for an angle of attack of 2° and a coning-rate parameter,' 4>£/V, of 
0.2. At this low angle of attack, good agreement is expected between the 

side-moment coefficient C n . and the planar damping-in-pitch coefficient, 

b 

as predicted by Eq. (2.39). Results of the computation are presented in 
Table 1 and are compared with the corresponding planar coefficients obtained 
from the unsteady inviscid flow- field results of Brong. 29 The agreement 
between the two methods is excellent with the maximum difference seen to be 
less than 2 percent. This demonstrates the capability of the present 
method to yield accurate results at hypersonic as well as at low supersonic 
Mach numbers. 


Table 1.- Comparison of Force and Moment Coefficients due to Coning with 

Corresponding Planar Coefficients; 10° Half-Angle Cone, M = 10.0, 

a = 2.0°, 4>£/V = 0.2, £ /£ = 0.6l. 

eg 



C Z 

/ \ 

C 

/ m \ 

/ r ^ iv 



(« c z>>) 

v a ' 

K w ) 


(*[V‘ J+Y Hv 

Coning 

-0.06717 

-0.00588 

0.01623 

-0.00626 

(Planar) 29 

-0.0671 

-0.00586 

0.0l6*+ 

-0.0063*+ 



3.6.3 Elliptic Cone at M - 2 

Computations vere carried out for a cone of elliptic cross section in 
coning motion at M = 1.97. The ratio of the major to minor axes of the 
cone was 3:2, and the base area was equal to that of a 7*75° half-angle 
circular cone of the same length. The angle of attack was fixed at for 
^ = 0 and tt/2, and the coning- rate parameter ranged from 0 to 0.15. Com- 
putations vere also carried out for the equivalent circular cone in coning 
motion. 

3.6. 3.1 Normal-Force Coefficient 

The normal-force coefficients obtained from the computations with the 
bodies in coning motion are presented in Fig. 18(a). Again, at fixed angle 
of attack no change in the normal- force coefficient was obtained as the 
coning-rate parameter was varied. Also shown in Fig. 18(a) are experimen- 
tal values of the normal-force coefficients of these bodies measured at 
zero coning rate by Jorgenson, 80 The good agreement between the computa- 
tional results and the experimental measurements demonstrates the capabil- 
ity of the method to compute the flow field about bodies of noncircular 
cross section. The method can be extended to bodies of arbitrary shape 
(i.e., airplane-like configurations). 

3.6. 3 . 2 Side-Moment Coefficient 

Computed values of the side-moment coefficient C n (°°;5,t) are P re “ 

$ 

sented in Fig. 18(b). Also shown is the theoretical value (Tobak and 
Wehrend 18 ) of the linear planar damping-in-pitch derivative 

(°°)+YC m .] for the equivalent circular cone. As before, for the 

^B S 

circular cone C n H 0 and there is good agreement between the computed 

$ 

value of C n and the linear damping-in-pitch derivative, as predicted by 
Eq. (2.39)* For the elliptic cone, C n does not vanish but can probably 
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be neglected in comparison with C n . . Under these conditions, Eq. (2.33) 

<t> 

indicates that C n . («; 6 ,ir/2) = 6[(L (°°) + YC m .] and 

<£ a 

C n (°°;<5,0) - 6 [ C n (°°) - yC n .]. These are the planar damping coefficients 

4 r B 3 

for oscillations in the directions indicated by the double-headed arrows 
in Fig, l8(b). The computational results confirm the expectation that 
oscillations parallel to the major axis of the elliptical cone should be 
less strongly damped than those parallel to the minor axis. 

4. CONCLUDING REMARKS 

Concepts from the theory of functionals have been applied to develop 
nonlinear formulations of the aerodynamic force and moment systems acting 
on bodies in arbitrary, slowly varying, large- amplitude motions. The anal- 
ysis, which proceeds formally once the functional dependence of the aero- 
dynamic reactions upon the motion variables has been established, ensures 
the inclusion of all pertinent aerodynamic terms within the resulting 
foimulations. The results confirm formally the intuitive notion that the 
instantaneous reactions on a body in an arbitrary motion can be compounded 
of the contributions from several characteristic motions: a steady motion 

with the instantaneous large values of the motion variables held fixed, 
and suitable perturbations about the steady motion. 

The analysis indicates that coning, where the nose of the body 
describes a circle around the velocity vector, characterizes the nonplanar 
nature of a general motion. With this motivation, a numerical method was 
developed for computing the inviscid flow field surrounding a body in con- 
ing motion. Computations carried out for circular cones in coning motion 
both at low supersonic and hypersonic Mach numbers confirm the adequacy of 
a linear formulation at low angles of attack. At larger angles of attack, 
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however, the forces and moments become nonlinear functions of the angle of 
attack* Computational results for the reactions on -the circular cone at 
the higher angles of attack agreed veil vith experimental measurements 
within the range of variables investigated* This indicates that the ini- 
tial nonlineai* behavior of the aerodynamic forces and moments is governed 
primarily by the inviscid flow and not, as previously postulated, by the 
viscosity-induced leeward-side vortices. 

Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 9^035 August 24 # 1973 
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APPENDIX A 


The formulation developed in the body- fixed axis system is 


p B f ' - 


C. (t) = C (“;d(t),g(t),V(t)) + C k (»ja(t),B(t),V(t)) 

K rL Pn 


B 


+ -^C> (« S a(t),B(t) f V(t)) + -§“ c k (»;a(t),g{t),V{t)) 

R q B R B 


* • 

+ |^e k (a(t) t B(t),V(t)) + |£c k .(a(tU(t),V(t)) 

R a R (3 

+ ““ C k -(a(t) ,£(t) ,V(t)) ; k= £,m,n (Al) 

t r ^ * 


while the formulation developed in the aerodynamic axis system is 


C k (t) = C k (“;6(t) ,ip(t) ,v(t) ) 


+ (“;$(t) ,i|t(t),V(t) ) 

R X 


+ ^C k (»;6(t),i|,(t),V(t)) + “• C k («;6(t) # *(t),V(t)) 

R q R 

+ |^C k U(t),*(t),V(t)) +1^^ (6(t)^(t),V(t)) 

R 8 R if; 

+ C k « ( 5 ( t ) ,i|/(t) ,V(t) ) ; k= £,m,n (A2) 

V R 

The moments expressed in the two axis systems are related to each other 
through 

C £ = (A3a) 

C + iC = e 1 ^(C + iC ) (A3b) 

m n m n 

while the motion variables in the two axis systems are related by 

P B = X + i|j (Alta) 
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q. B + ir B = e -1 ^(q + ir) 


(kkb) 

a + ig = 6 e llfl (A 4 c) 

Substituting the results of Eq. (a 4 ) into Eq. (Al), we obtain 

C,(t) = C. (-;a(t),|(t),V(t)) (»;a(t) ,g(t) ,V(t)) 

* V R P B 

+ v^l^k (*;o(t) ,p(t),V(t))cos 4 ; - C k (®;a(t ) ,g(t) ,V(t) )sin 4*1 

RL q B r B J 

+ (°°;a(t) , 0 (t) ,V(t) )sin \p + C k (® ;a(t ) ,g (t) ,V(t ) )cos 4 * I 

RL * q B r B J 

+ l^fSk. (a(t) , 6 (t ) ,v(t) )cos 4 * + C^. (a(t) ,£(t) ,V(t))sin 4*1 
RL a § J 

+ 6 % 1 -Ofc. (a(t) ,g(t) ,V(t) )sin \p + C k , (a(t) ,§(t) ,V(t) )cos 4 * 

La & 

+ j C k p (°°;a(t) ,$(t) ,V(t))l + ~ C k ^.(a(t),g(t),V(t)) 

B -* V_ 


; k = £,m,n 


(A5) 


When the results of Eq. (A5) are substituted into Eq. (A3), the following 
equivalences between coefficients formulated in the two-axis systems are 
obtained: 


X 

4 ( 


■ B 


d B 


i B 


6 a 


A A ^ 

■ja.B.V) 

(A6a) 

1 

(A6b) 

(«°;a ,3 ,V)cos ip — (®;a,£,V)sin tp 

i r B 

(A6c) 

(®;a,g,V)sin 4 + {»;a,f3,V)cos 4> 

i r B 

(A6d) 

aj6jV)cos i|) + (a,|,V)sin ip 

6 

(A6e) 

a,B,V) 

(A6f) 
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and 


Cjj, (6,4> S V) = Cjj, (»;a,3,V) + 5[-C£. (a,B ,V)sin 4* + C^, (a,B ,V)cos 
p b S. B 

(A6g) 


c m («j5,4/ t v) 

X 

Cn^(-;«,v»,v) 

Cm r (“;6,t s V) 

5 

Prny('S»4 , > v ) 

^m. ( 6,4 l s^) 
4 


C m (“;a,3 1 V)c°s 4> - C^(«>;a s B,V)sin ^ 


Cjj^ } V)cos (®;a,B,V)sin 4> 

C m («;a,3,V)cos ip + C n («;a,B ,V)sin 2 4; 
q B r B 

— [Cjn (»;a, 8 ,V) + C n (°° ; a ,3 ,V) ]cos 4) sin ip 
r B q B 

Cm (®;a,B,V)cos 2 4> - C n («>-,a,B>V}sin 2 i|j 
^B q B 

+ [Cm (°°;a,B,V) - C n (°°;a , 6 ,v) jcos 4» sin <p 
q B r B 

C m . (a>3 ,V)cos 2 4» - fi n . (a,B,V)sin 2 4j 
a B 

• ( ct 9 B ) — Cj-j . ( ct 3 $ jV ) ]cos 4 » sin 

^ A 

3 a 

C m ^(a,8,V)cos jJj - C n ^.(a,B,Y)sin 4/ 

Cj,^ (°°;a,B,V)cos ip - («>;a,B ,V)sin i|/ 

+| 5{C m . (a»B »V)cos 2 4j + C n ^(a,B»V)sin 2 ^ 

3 a 

-[C m .(a s B,V) + C n . (a,B,V) ]sin ^ COS 
a 3 


(A7a) 

(ATb) 

(A7e) 

(A7d) 

(A7e) 

(A7f) 


(A7g) 


The correspondences for the coefficients in C n (t) may he obtained from 

/ \ A A A A 

Eq. (AT) by replacing with and C with ~C^ # Similarly, the anal- 
ogous matches for the force coefficients may be obtained by replacing 
£,m,n with X,Y,Zj respectively, in Eqs. (A 6) and (A7)* 
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s a(£) = a(r) + a(r)(£-r) 


^q B (£) = q 0 (r) + q B (r)(^-r) 
















^^tip sin <r 



(a) Coning motion. (b) Conical flow approximation 


Figure 10.- Approximation of initial solution. 
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Figure 12*- Concluded 
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Figure 17.- Typical vapor screen photographs; 10° half-angle cone, 

Mach number = 2.0. 
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(a) Normal -force coefficient. 


Figure l8.- Normal-force and side-moment coefficie 
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